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Abstract. The stable propagation of jets in FRII sources is remarkable if one takes into account that large-scale jets are 
subjected to potentially highly disruptive three-dimensional (3D) Kelvin-Helmholtz instabilities. Numerical simulations can 
address this problem and help clarify the causes of this remarkable stability. Following previous studies of the stability of 
relativistic flows in two dimensions (2D), it is our aim to test and extend the conclusions of such works to three dimensions. 
We present numerical simulations for the study of the stability properties of 3D, sheared, relativistic flows. This work uses a 
fully parallelized code (Ratpenat) that solves equations of relativistic hydrodynamics in 3D. The results of the present simula- 
tions confirm those in 2D. We conclude that the growth of resonant modes in sheared relativistic flows could be important in 
explaining the long-term collimation of extragalactic jets. 
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1. Introduction 

Extragalactic jets from AGN present a r norphological di- 
choto my between FRI and FRII type jets ( Fanaroff & Rilev , 
\l974 ). The former (e.g., 3C 31. lLaing & Bridlell2002l) show 
a disrupted st ructure at kiloparsec sc ales, whereas the latter 



(e.g., Cyg A, ICarilli & BarthelL 1 19961) are highly collimated. 



Although the origin of this dichotomy is a complex combi- 
nation of several intrinsic (e.g., jet power) and external (i.e., 
environmental) factors, the stable propagation of jets in FRII 
sources is remarkable considering that large-scale jets are sub- 
jected to potentially highly disruptive three-dimensional (3D) 
modes of the Kelvin-Helmholtz (KH) instability. Any pertur- 
bation may couple to an instabili ty and grow in am plitude, be- 



erties of the KH instability in relativistic cylin drical jets. Th e 
effects of a sh ear layer have been examined bv lUrpinl (120021) . 
More recentlv. 'Harded (12008) has studied the stability of mag- 
netized relativistic jets with poloidal magnetic fields. 

An extension of the linear stability analysis in the relativis- 
tic case to the weakly nonlinear regime has been performed by 
Hanaszl (11995. .1997.) and led to the conclusion that KH insta- 



bility saturates at finite amplitudes because of various nonlinear 
effects. The most significant effect results from the relativistic 
character of the jet flow, namely from the velocity perturba- 
tion not exceeding the speed of light in the reference frame 
of the jet. This result was confirmed by numerical simulations 
dPerucho et al.li2004a_b.) . 

In the purely non linear regim e, Marti et al.l (Il997h 



X 



importance of studying their development and growth from the 

linear to the nonlinear regime in terms of the flow parameters. 

Jet stability in the relativistic regime has been thor 

oughly s tudied in different scenarios and parameter ranges (see 



> coming potentially disruptive (Perucho et all |2005|) , hence the JHardee et al .[ (M and lRosenetal.1 (M demonstrate thai 



Hardeel 12006 , for a recent review on the state-of-art), from 



both the analytical and numerical perspectives. The linear anal- 
ys is of KH instability i n relat ivist ic flows started with the wor k 
of iTurland & Scheueil 11976") and'Blandfo rd & Pringld (1 19761) . 
who derived and solved a dispersion relation for a single plane 
boun dary between a relat ivisti c flow and the a mbient medium. 
Next, Ferrari et alj (119781) and iHarded (119791) examined prop- 
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jets with high Lorentz factors and high internal energy are in- 
fluenced very weakly by the Kelvin-Helmholtz instability, con- 
trary to the cases with lower Lorentz factors and lower internal 
energies. The former do not develop modes of KH instability 
predicted by the linear theory, which is interpreted as the result 
of a lack of appropriate perturbations (triggered by the back- 
flow in this kind of simulations) generating the instability in 
the system, because in the limit of high internal energies of the 
jet matter the Kelvin-Helmholtz instability is expected to de- 
velop with the highest growth rate. 

A combination of linear and nonlinear studies of jet sta- 
bility in the relativistic case was applied for the first time by 
Hardee et al.l (11998) in the case of axisymmetric, cylindrical 
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jets and then extended to the 3 D case by iHardee et al. (1200 Ih 
in the spatial approach and by IPerucho et alj (l2004allbl) in the 
temporal approach, for 2D sl ab and cilyndrical jets , including 
the case of sheared flows in IPerucho et alJ OOOSh . Recently, 



Mizuno et alJ (l2007h have extended this study to 3D relativis- 



tic magnetized jets, which may be stabilized (with respect to 
current-driven modes) by a sheath layer surrounding the jet 
core. 

Focusing on the com bination of linear and nonlinear stud- 



Perucho et al.l (12004a b) studied the growth of single sym- 



les, 

metric modes in relativistic jets, sweeping diff'erent parame- 
ters, such as jet temperature, Lorentz factor, and density ra- 
tio, with the external medium. The linear work confirmed that 
KH modes grow faster in slower and hotter flows than in cold 
and fast ones. In Perucho, Marti & Hanasz (2005, Paper I from 
now on), this work was extended to sheared jets to study the 
growth of a mixture of symmetric and antisymmetric modes in 
this case. In IPerucho et al.l (|2004b) and Paper I, the authors de- 
termined the stability properties of relativistic flows for a wide 
range of parameters in density contrast (lO"'* - 0.1), Lorentz 
factor (5 - 20) and specific internal energy (0.08 - 76.5 c^). It 
was shown that the inertia of the jet, in terms of the density 
contrast with the external medium, the relativistic Mach num- 
ber and the Lorentz factor are the parameters that determine the 
stability of jets in the nonlinear regime, as a relativistic coun- 
terpart to the work of iBodo et al.l (119941) in th e classical case. 
Moreover, in Paper I and Peruch o et al.l (l2007b (Paper II, from 
now on), it was shown that the growth of high-order Kelvin- 
Helmholtz modes developing in the shearing layer, hereafter 
referred to as resonant modes, could dramatically change the 
nonlinear evolution of the flow in the case of cold and fast jets. 

The results could be summarized as follows, a) Cold and 
slow jets (small relativistic Mach number and Lorentz factor) 
are unstable with the growth of the instability and are disrupted 
after the generation of a shock front crossing the boundary be- 
tween the jet and the ambient medium, which results from the 
development of long wavelength KH instability modes (USTl 
jets), b) Hot and slow jets (intermediate values of relativistic 
Mach numbers and Lorentz factor) are also unstable, but dis- 
rupted and mixed in a continuous way by the growth of the mix- 
ing layer down to the jet axis (UST2 jets), c) Faster (high val- 
ues of Mach number and/or Lorentz factor) jets develop short- 
wavelength, high-order modes that grow in the shear layer and 
saturate the growth of the instability without loss of either col- 
limation or mixing, and generate a hot shear layer around the 
core of the jet (ST jets), d) Lighter jets are more unstable than 
denser ones (USTl). 

In this work, we extend the previous studies to 3D and show 
the existence of a new mechanism to prevent jet decoUimation 
based on the development of resonant KH perturbations. With 
the aforementioned experience in mind, and taking into account 
that 3D simulations are highly demanding in terms of compu- 
tational resources, we present here three numerical simulations 
using parameters that are representative of the different stability 
regions found in Paper I, but now including helical and elliptic 
modes. It is our aim to test whether the same conclusions are 
valid and, thus, if resonant modes provide a general stabihsing 
mechanism in the absence of strong magnetic fields. 



The paper is organized as follows. Section |2] is devoted to 
presenting of a new 3D-RHD code for these simulations and the 
numerical setup used. Section |3]contains the results of the lin- 
ear analysis for the parameters considered and the simulations. 
In Section |4] the results are discussed and the conclusions of 
this work are presented in Section |5] 

2. Numerical simulations 

This work was done with Ratpenat, a new finite-difference code 
based on a high-resolution shock-capturing scheme that solves 
the equations of relativistic hydrodynamics in 3D written in 
conservation form. It is based on the 2D code used in all the se- 
ries of works on jet stability performed in our group (see Paper I 
and refer ences t herein), and shares algorithm and structure with 
Genesis (lAlov e t al., 1999). Ratpenat has been parallelized us- 
ing a hybrid scheme with both parallel processes (Message 
Passing Interface, MPI) and parallel threads (OpenMP) inside 
each process. This hybrid scheme is well adapted to the typical 
supercomputer architecture, i.e., a cluster of multi-processors. 
Ratpenat can exploit all the memory available in a compu- 
tational node and keep all the cores in the node busy with 
OpenMP threads. 

The parallelism between computational nodes is exploited 
with the MPI layer. Moreover, this hybrid scheme reduces the 
total amount of communications, because the number of MPI 
processes is reduced to the number of nodes and prepares the 
code for future supercomputer environments where load imbal- 
ance between MPI processes can be solved dynamically with 
the proper management of threads priorities. The distribution 
of computing load among the different participating processors 
is done by dividing the numerical grid in domains along the 
jet axis and distributing each portion to one of those proces- 
sors. Communication is made before every temporal subcycle 
and the data provided are used as boundary conditions for each 
processor. The super-computer Mare Nostrum, at the Barcelona 
Supercomputing Centre (BSC), is composed by nodes that in- 
clude four processors with shared memory. The code was tested 
to scale up to 256 nodes in Mare Nostrum, i.e., 1024 proces- 
sors, with an appropriate grid size (xxyxz = 256x8192x256). 

The initial conditions in the simulations are similar to those 
used in Paper I for the 2D case: a portion of an inifinite jet, 
represented by periodic boundary conditions in the axial direc- 
tion and outflow boundary conditions in the radial directions. 
The size of the grid is 16 Rj (jet radii) along the axis and 8 Rj 
transversally (4/?j on each side of the jet axis), plus an ex- 
tended grid up to 200 Rj on each side, resulting in a grid with 
512 X 512 X 512 cells. The resolution is 32 cells/Rj in the 
homogeneous grid. This is much less than that in the 2D sim- 
ulations. The influence of the resolution on the results will be 
discussed in the next section. 

Each jet is initially in pressure equilibrium with the ambient 
medium. A shear layer is introduced between the jet and the 
ambient medium for axial velocity and rest-mass density: 



vAr) = 



Po(r) = POa - 



cosh(r'") 

POa - POj 



cosh(r™) ' 



(1) 



(2) 
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where r is the radial coordinate, poj and poa are the density 
of the jet and the ambient medium, respectively, and m is the 
steepness parameter of the shear layer In the 3D simulations, 
we used m - 16, which implies a wider (~ 0.3 Rj) layer than in 
the 2D jets in Paper I (~ 0.2 Rj). Both the jet and the ambient 
medium are described as relativistic ideal gases with adiabatic 
index F = 4/3. Thus, we implicitly assume that the jet is sur- 
rounded by a hot cocoon. 

We chose three representative models out of those studied 
in Papers I and II: Model BOS, DIO, and B20, belonging to 
each of the three families mentioned in the Introduction. Their 



parame ters are summarized in Table |2] (see also IPerucho et al 



(l2004ah and Paper I). The columns give, for each model, the 
Lorentz factor, the relativistic Mach number, the density con- 
trast, the adiabatic index, the specific energy density of the jet 
and the ambient medium and the sound speed of the jet. Models 
BOS and B20 are shown along the paper as extreme cases, be- 
longing to groups UST 1 and ST in the stability plane derived 
in Paper I (Fig. 21 in Paper I), so these models are taken as 
paradigms of those two groups. Model DIO shares properties of 
both BOS and B20, belonging to the UST2 group in the stabil- 
ity plane. Following the similarities for all models within these 
groups, we expect that any conclusion derived for the mod- 
els presented here can be generalized to all the models falling 
in the same regions, including lighter, slower, faster, hotter, or 
colder jets depending on the group. 

A combination of linear pinching, hehcal, and elliptic sinu- 
soidal perturbations is excited in all three simulations. These 
represent the fundamental mode of the numerical box {A = 
16 Rj) and the first three of its harmonics. The amplitude of 
the pressure perturbation is set by taking into account that the 
linear regimes are longer for colder and faster jets and that the 
initial amplitude does not change the long-term evolution of 
the jet (see Paper I). The amplitude of the perturbation is thus 
greater in model B20 in order to avoid a very long linear regime 
(see, e.g. iPerucho et al.l l200Sb . In Paper I it was proven that 
these perturbations couple to the fastest growing modes at the 
given wavelength and that any other shorter harmonics, such as 
the short, resonant modes, can be excited. Therefore, we have a 
representative sample of the different types of modes that grow 
in such a scenario and are able to test the general response of 
each set of parameters to the growth of the instability. 

The simulations were run for between 18 and 30 days in 32 
nodes (128 processors) in Mare Nostrum (BSC), depending on 
the run. The number of nodes used was chosen on the basis of 
the grid dimensions, as explained above. 



3. Stability of three-dimensional relativistic jets 

3.1. Linear regime 

The solutions to the linear problem are shown in Fig.[T] These 
solutions (for the wave-number and frequency of an instability) 
were obtained for the differential equation resulting from the 
linearization of the equations of relativistic hydrody namics for 
a sheared jet in cylindrical coordinates (e.g.. lHardeei. 2000. for 



the vortex sheet case, and Paper II for a sheared, slab geometry 
jet): 

dp dp 

^+Fi{r)^+F2{r)p^Q, (3) 

or-^ or 

with r the radial coordinate, p the perturbed pressure, 

1 dv^fi ylvzd'^ - kzVzfi) + kz dpydr 

Fi(r) = 2 ; , (4) 

r dr CO - k,Vzfl Po^o 

and 

Fiir) = — ^--+(r^z.o(w - k,v,^ + k,) {u>v,fl-k-X (S) 

where y is the Lorentz factor, v, the axial velocity, u the com- 
plex frequency, k, the complex wavenumber, n the azimuthal 
wavenumber (0 for pinching modes, 1 for helical modes and 
2 for elliptic modes), p' the proper energy density, h the spe- 
cific enthalpy, c^j the sound speed of the jet and the subscript 
refers to the equilibrium values. The radial derivatives take 
a shear layer into account in the axial velocity and the energy 
density. 

Equation (O was solved using the shooting method. This 
implies integration with a Runge-Kutta of variable step, from 
the axis to a given radius far from the jet boundary, for each pair 
of temptative values for w and k~. The boundary conditions on 
the jet axis are fixed in terms of the symmetry properties of 
the perturbation. The roots are found by fixing a real value of 
the wavenumber and looking for the roots (complex frequency) 
of the boundary condition that states that the amplitude of the 
wave must be zero far from the jet. The root-finder used with 
this purpose is the Miiller method for equations of complex 
variable (check Paper II for details). This process of calcula- 
tion explains the point-to-point structure of the results shown. 
In addition, Bessel functions with complex arguments, which 
appear as solutions to eq. (O, are numerically computed. This 
is the first time that these kinds of results are systematically 
obtained for relativistic, sheared, cylindrical jets. 

The results shown in Fig. [1] reveal that the helical and el- 
liptical surface modes dominate the growth rates at the longer 
wavelengths (low-order modes) for the three models. The plots 
show that the surface mode is damped at short wavelengths 
for colder jets (BOS and B20, top and bottom panels) when 
compared with the hotter jet (DIO, central panels). The growth 
rates tend to decrease from BOS and DIO to B20, i.e., with 
the Lorentz factor and relativistic Mach number. The resonant 
modes dominate the solutions at the shorter wavelengths in all 
the models. Their growth rates are higher than those of the 
low-order modes. These results coincide with those obtained 
in Paper II for jets with slab geometry. 

We excited several modes for each model {k ^ 
0.4, 0.8, 1.6, and 3.2 /Jji). In the case of BOS, the pinching 
perturbations excite wavelengths surrounding the maximum 
growth rate of the first body mode, of the second body mode, 
and the resonant modes, whereas the helical and elliptical per- 
turbations excite the surface, first body, second body, and reso- 
nant modes. In the case of DIO the pinching perturbations trig- 
ger the surface (at a very low growth rate), first and fourth body 
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Jj 


M, 


PjlPa 


r 


£/c2) 


Sa{c^) 


Csjic) 


resolution (cells/i?^) 


B05 (3D) 


5.0 


13.2 


0.1 


4/3 


0.42 


0.042 
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32 / 64" 


BOS (2D) 


- 


- 


- 


- 


- 


- 


- 


32/128x256* 


DIO (3D) 


10.0 


14.2 


0.1 


4/3 


60.0 


6.0 


0.57 


32 


B20 (3D) 


20.0 


54.0 


0.1 


4/3 


0.42 


0.042 


0.35 


32 


B20 (2D) 


- 


- 


- 


- 


- 


- 


- 


32 / 128x256* 



Table 1. Parameters of the simulations discussed in this paper. " This simulation was repeated (in the linear regime only) doubling 
the resolution. * The numbers correspond to the axial times radial resolution when it is not the same in all the dimensions. 
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Fig. 1. Solutions to the linear problem. Upper panels show the solutions for model BOS, central panels show the solutions for 
model DIO, and the bottom panels show the solutions for model B20. The left column displays the pinching modes, the central 
column shows the helical modes, and the right column shows the elliptic modes. In each of the panels, the upper lines represent 
the real part of the frequency, whereas the lower lines represent the imaginary part of the frequency, i.e., the growth rate. The 
surface modes are those with no zeros between the jet axis and its surface, whereas the body modes have one or more zeros in 
this region. They appear from left to right in the images. 



modes, and the helical and elliptical perturbations excite the 
surface and fourth and third body modes, respectively. Finally, 
in the case of B20, all the perturbations excite the first or sec- 
ond body modes and the resonant modes, as the surface mode 
is damped at the low values of the wave-number {k < QAR.^ 
for the helical surface mode and k ^ 0.4 /?t' for the elliptical 
surface mode). We can thus anticipate from these results that 



the resonant modes may influence the growth of perturbations 
in all three cases, but will be more important in the case of 
pinching and helical modes in B20, due to their higher relative 
growth rates. 

During the linear evolution, the resonant modes are excited 
either directly or as harmonics of the longer wavelengths and 
may grow faster than the latter, as explained in the previous 
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Fig. 2. Left panels: Pressure perturbation across the jet axis for models BOS (left panels) and B20 (right panels) at two different 
times (f I and f2 in Fig. [3] for BOS and f i and t^ for B20), within the linear regime (upper row) and close to the end of the linear 
regime (lower row). Arrows in Fig. [3] indie ate these times. Right panels: Specific internal energy across the jet axis for models 
BOS (left panels) and B20 (right panels) at the same times. The two lines (solid and dotted) stand for two different positions along 
the jet axis, at half the grid-size distance to each other, where the cuts are done. Note the change in scale between the upper and 
lower plots. 



paragraph. In Paper II, the authors showed that these modes 
grow faster in the shear layer between the jet and the ambient. 
Thus if the resonant modes dominate the growth, the largest 
mode amplitudes are expected towards the shear layer. The up- 
per panels in Fig. |2] show the typical profiles produced by the 
growth of resonant modes (see Figs. 4 and 9 in Paper II) for 
models BOS (left) and B20 (right). The radial profile of the 
perturbations in these jets is clearly dominated by the reso- 
nant modes. However, this does not imply that they have the 
highest amplitude in the whole volume of the jet, as revealed 
by a Fourier analysis along the jet (not shown here), which 
shows that they are in competition with the longer wavelengths, 
mainly in the case of BOS. Axial profiles of the pressure pertur- 
bation (not shown) taken at r = 0.5 Rj reveal long antisym- 
metric structures along the jet for model BOS, together with 
the short scale oscillations, whereas in the case of B20, the 
long structures appear to be mainly symmetric (due to the high 
growth rate of the elliptical modes) and short scale oscillations 
of comparable amplitude are also visible. 



The lower panels in Fig. |2] show the same variables at a 
later time, close to th e end of the linear regime (as defined in 
Perucho et al.Ll2004 aV The evolution of the pressure perturba- 
tion in time is very similar, in both cases, to that shown in Fig. 
9 of Paper II, with the larger amplitudes initially showing up in 
the sheared region. We show, in the right panels of Fig. [31 two 
cuts of specific internal energy across the jet for both models, at 
the same times as in the left panels. Both the upper (in the lin- 
ear regime) and lower (about the end of the linear regime) plots 
reveal that the resonant modes convert more energy into heat 
in the shear layer in the case of model B20. Thus, as expected 
from the linear solutions (see Fig.[TJ, these modes influence the 
dynamics of this jet more. 



3.2. Nonlinear regime 

In previous works the stability of jets in the nonlinear regime 
was characterized in terms of mixi ng and transfer of linear mo- 
mentum to the ambient medium (IPerucho et al.l l2004bl and 
Paper I). Following these criteria, we discuss here the stabil- 
ity of the simulated jets. The axial momentum in the simulated 
jets is displayed versus time in Fig.[3j normalized to their value 
at f = 0. The longitudinal momentum in the jet is computed by 
multiplying the actual axial momentum in each cell by the jet 
mass fraction. During the linear regime, we observe no transfer 
of linear momentum from the jet to the ambient. This regime 
is shorter in these simulations than in the 2D ones because the 
initial amplitudes of the perturbations have been set to be larger 
with the aim of reducing the computational times. In the non- 
linear regime, the plots show a clear trend following the results 
reported in Paper I: model BOS loses the largest amount of axial 
momentum and model B20 keeps a large part of it at the end of 
the simulation. The latter was kept running longer than mod- 
els BOS and DIO to check that this trend is valid during long 
times in the nonlinear regime. Model DIO was stopped at an 
earlier stage than BOS because the time-step became progres- 
sively shorter during the nonlinear regime and the simulation 
demanded too much time to be continued from this point. 

Figures m (model B0S),|5](D10), and|6](B20) display axial 
cuts at the final snapshot in Lorentz factor and tracer (jet mass 
fraction) for the three simulations. The maximum Lorentz fac- 
tor in the maps of models BOS is y = 1 .67, in the maps of model 
DIO is r = 13.1 - 13.6, and y = 21.1 - 21.7 for model B20. 
Regarding the tracer, white indicates original, unmixed jet ma- 
terial. From these maps, we see that BOS is completely mixed 
and DIO is strongly entrained, although both models keep a 
central region with the higher values of this variable, whereas 
model B20 keeps a central spine with little or no mixing. 

The maps show the effect of the different modes in the 
structure of the flow and a gradient in the stability of the 
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jets from BOS and DIO to B20, as could be expected from 
Fig. |3] Although the mixture of growing modes makes it dif- 
ficult to identify the individual ones, these images allow us to 
see which wavelenghts dominate the nonlinear regime. Model 
BOS shows a helical structure with the size of the box that can 
be attached to the surface mode (see Fig. [TJ, but also small- 
scale efficient mixing structures resulting from the growth of 
the short-wavelength resonant modes. The Fourier transform 
shows competition among different harmonics, mainly the fun- 
damental mode of the box, followed in amplitude by the first 
harmonic and a mixture of higher order harmonics, which co- 
incides with the predictions from the linear regime. The rapid 
disruption of this jet is the consequence of the distortion of the 
jet surface by the longer wavelengths. 
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Fig. 3. Evolution of axial momentum in the jets with time. The 
solid line stands for model BOS, the dotted line for model DIO 
and the dashed line for model B20. Indicated with arrows, the 
times at which the cuts in Fig. |2] were taken: fi is the time of 
first cut for both models, f2 the time of the second cut for BOS, 
and f3 the time of the second cut for B20. 



The central spine of model DIO shows a short wavelength 
(1 - 2Rj), helical structure, which could be identified with the 
fourth body mode (see Fig. [T]i, and, on the same scales, expan- 
sions and contractions of the jet, which could come either from 
the fourth pinching body mode or from the third elliptical body 
mode. These structures are surrounded by a layer where small- 
scale mixing occurs. In this region, an 8 Rj helical pattern can 
be seen. This result clearly shows that different wavelengths 
may show up at differ ent jet radii, generating a radial profile, 
as already shown in Perucho et al.l (12006) . The Fourier trans- 
form indicates that the highest amplitudes correspond to the 
fundamental mode and first harmonic at the end of the linear 
regime, followed by high-order harmonics, which keep grow- 
ing slowly after the end of the linear regime until saturation of 
the growth of the radial velocity and domination of the nonlin- 
ear regime in the long term. The growth rates shown in Fig. [1] 
agrees with the importance of the third to sixth body modes 
during the evolution. However, in contrast to the colder mod- 
els, the linear solution does not show small-scale resonances 



for high-order body modes, but the resonances occur for those 
low-order body modes, which may be the reason for the final 
disruption of the jet. The structure of jet DIO, with maxima 
and minima in the Lorentz factor along the axial directions will 
make shock waves arise naturally, which may mix the jet fur- 
ther and eventually disrupt it. 

Model B20 presents a more collimated structure and shows 
a combination of expansions and contractions plus a helical 
pattern with the wavelength of the box, corresponding to the 
first body helical mode. As in model DIO, there is a thick mix- 
ing layer surrounding the jet core. The Fourier transform also 
shows that the fundamental mode of the box (corresponding to 
the helical first body mode) dominates during the linear regime, 
but also the fast growth of the higher harmonics, which, like 
before, surpass the amplitude of the long wavelength during 
the time interval between the end of the linear regime and the 
saturation of the radial velocity. This competition between the 
longest wavelength and the shorter ones can be deduced from 
Fig. [U In the case of model B20, the structure at the end of the 
simulation points towards further mixing, although this jet has 
a higher Lorentz factor and a more uniform structure than DIO 
at the end of the simulation. 

Figure|7]shows transversal cuts at halfway through the grid 
at the end of the simulation. The figure shows the final dis- 
torted structure and turbulent mixing in all three models. The 
top panels show that model BOS is completely mixed (see the 
tracer in the right panel), whereas DIO and B20 still keep a 
fast, unmixed spine, which is larger in the latter. B20 shows a 
hot shear layer including the slower, outer regions of this un- 
mixed core, contrary to DIO and BOS, where the temperature 
decreases outwards. The centroid of model B20 is displaced 
left and down from the centre of coordinates, indicating helical 
motion in competition with the elliptical modes, responsible for 
the oblique elongation of the jet section. There are hints of a 
similar deformation in model DIO, but no significant displace- 
ment of the centre of the jet is detected in this case. This may 
be because the dominating helical structure is given by a short 
wavelength mo de, which only generates small displacements 
(lHardeel l2000). Also model BOS shows an elongation (vertical 
in this cut), but any hints of regular structure have been smeared 
out by the mixing and disruption of the original flow. It is also 
difficult to assess whether the axis of the jet is displaced by 
helical modes; this can be better observed in the axial cuts. 

To have a global view of the parameters along the jet axis, 
we show in Fig. [8] mean axial values of different physical pa- 
rameters for the three models at the end of each simulation. 
In these panels, the displacement of the axis due to the helical 
modes is promediated and thus we always find the centre of the 
jets around their original position. The left panels display the 
Lorentz factor and show barely relativistic motion for model 
BOS, with maximum Lorentz factor y - 1.36, as compared to 
models DIO (y = 8.28) and B20 (y = 18.4). The central panels 
show the logarithm of the specific internal energy and show the 
same behaviour as observed in the 2D simulations of Papers I 
and II: the generation of a hot shear layer around the core in 
those models in which the resonant modes already dominate 
during the linear regime (model B20), due to the strong dis- 
sipation of kinetic into internal energy in this area, when the 
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Fig. 4. Axial cut of Lorentz factor (left) and jet mass fraction (tracer, right), along the central planes for model BOS at the end of 
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Fig. 5. Axial cut of Lorentz factor (left) and jet mass fraction (tracer, right), along the central planes for model D 1 at the end of 
the simulations. 



resonant modes achieve the nonlinear regime. The right pan- 
els show the tracer, which represents the percentage of initial 
jet material in a cell. In the case of model BOS, the whole jet 
is mixed with the ambient medium, as indicated by the maxi- 
mum original jet mass fraction in a cell (~21 %), model DIO 
shows some mixing down to the axis (with maximum mean jet 
mass fraction ~ 79% on the jet axis), whereas in model B20 the 
central spine remains unmixed and coUimated. Regarding the 
width of the mixing layer, we observe that model BOS presents 
the wider mixing region, whereas DIO has the thinnest one. 
However, if we consider that DIO was stopped much before 



the other models and that, from Fig. [3] the expected evolution 
is close to that of BOS, we expect that the mixing layer of DIO 
after some time should be wider than that of B20. 



4. Discussion 

4.1. Effects of the numerical resolution 



In lPerucho et al.l (l2004ah it was shown that the resolution used 
in the simulations is crucial to reproducing the linear growth 
rates of the modes in the linear regime. In particular, the resolu- 
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the simulations. 
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tion may affect the evolution of the shortest (resonant) modes, els BOS and B20 were performed in which we used the same 

which can be damped by numerical viscosity, and also have resolution as in the 3D simulations, and compared them with 

implications in the nonlinear regime with respect to, e.g., mix- the results obtained for cylindrical jets in 2D (see Appendix B 

ing. Thus, we performed several tests to check the influence of in Paper I). In these low-resolution simulations, we only ex- 

the resolution on our results. First, 2D numerical tests of mod- cited the pinching modes. In the case of BOS, the disruption 
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of the 2D jet is faster, both in the low-resolution and the high- 
resolution cases, than in the 3D case, because long-wavelength 
modes have the largest amplitudes at the end of the linear 
regime. In 3D, the growth of short-wavelength modes (heli- 
cal and/or elliptic) may be responsible for the slower mixing 
process. Between the 2D simulations, the low-resolution one 
shows a slower mixing and momentum transfer to the ambient, 
implying that the resolution affects this process, as expected. 
In model B20, we identify the effects of resonant modes also in 
the low-resolution 2D simulation, which result in a slow trans- 
fer of axial momentum from the jet to the environment during 
more than IQQQRjjc code time units. However, the amplitude 
of these resonant modes at the end of the linear regime is less 
than in the high-resolution case. This translates into a decrease 
in the Lorentz factor down to the jet axis and more transfer of 
momentum from the jet to the ambient in the low-resolution 
case, around an equivalent time to what is reached by the high- 
resolution simulation. We continued the low-resolution simu- 
lation further and found that the oscillations in the amplitude 
of all the (excited) m odes close to their limit of growth (see 
Perucho et all 12004a. and Paper I) end up in a long wavelength 
disrupting the flow. This occurs at time t ^ 2500 Rjjc, which, 
for a flow with radius, e.g., Rj - 100 pc propagating close 
to the speed of light, implies covering a distance greater than 
200 kpc before disruption. We should expect that higher resolu- 
tion would allow for an undamped growth of resonant modes, 
thus making these jets even more stable. 



In addition, we performed a 3D simulation of model BOS 
with twice the resolution used here. In this case, we would ex- 
pect a faster growth of the resonant modes already detected in 
the 3D low-resolution simulation. This simulation was run in 
the linear and postlinear regime only during ^ 150 Rj/c code 
time units due to the computational time limitations. The re- 
sults confirm the previous statements: the (helical or elliptic) 
resonant modes grow faster and make the mixing of this jet 
slower (up to the moment in which we stopped the simulation) 
when higher resolution is used. 

In conclusion, the low resolution may quantitatively affect 
the results, damping the growth of the resonant modes and also 
reducing the mixing and momentum transfer rates, but the re- 
sults presented here are qualitatively equivalent to those ob- 
tained in 2D with higher resolution. 

4.2. Implications for jet collimation 

The results are thus consistent with those of Papers I and II, 
where the authors showed that the growth of resonant modes 
in sheared jets can influence the evolution of relativistic flows, 
by generating small-scale mixing in the shear layer and avoid- 
ing the generation of large scale nonlinear structures that may 
disrupt the whole jet. This result could be important for our 
understanding of the properties of stability of relativistic out- 
flows on the kiloparsec scales. Even though the Lorentz fac- 
tors of both FRI and FRII jets appear to be similar on the par- 
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sec scales JGiovannini et aLLbOOl VCelotti & Ghisellinil 2008h 
the pr esent paradigm of FRI jets (Bicknell, Il984i iLaingl Il993 



1996) states that, on kiloparsec scales, these jets are deceler- 
ated by entrainment of gas and their properties become differ- 
ent from those of FRll jets on these scales. Two main pro- 
cesses have been invoked to explain the entrainment in jets 
and studied in the literature; (i) entrainment through mixing 
in a turbulent shear layer between the jet and the medium 
dPe YoungL Il986i 11993 : Bicknell, 'l99 4, and (ii) entrainment 
from stellar mass losses (.Komissarovi Il994i Bowman et all 
1996t iLaing & Bridlel l2002h . Recent numerical simulations 
deal with this last process (Perucho et al., in preparation). 
Concernin g the process of entrain ment through a turbulent 
shear laver. lPerucho & Martil (l2007h showed, via a 2D axisym- 
metric simulation, that a recollimation shock in a light jet, in 
reaction to a steep interstellar density and pressure gradients, 
may trigger nonlinear perturbations that lea d to jet mixing and 
disruption with the external medium. Later, iRossi et al.l (12008') 
obtained a similar result in 3D for a homogeneous ambient 
medium, where the shock is ind uced by an overpressured co- 
coon. Also M eliani et alJ (120081) studied the influence of a dis- 
continuous jump in the external density gas on the jet disrup- 
tion in hybrid FRl/FRll sources, where the jet and counterjet 
present different morphology, and this is adscribed to such en- 
vironmental irregularities. 

Regardless of the origin of the deceleration, the Lorentz 
factors in FRI jets after 1-2 kpc may already be small enough 
to favor the growth of KH modes, triggering the mixing pro- 
cess. In the context of the FRI/FRII dichotomy, the evolution 
of the KH instability can thus be seen as a consequence of 
this dichotomy and a cause of the kiloparsec-scale morphol- 
ogy of the two types of jets, where in one case the growth of 
the KH modes favor mixing and in the other may favor fur- 
ther collimation. In the case of FRII jets the strong recolli- 
mation shocks or strong entrainment discussed in the previ- 
ous paragraph do not seem to happen. Nevertheless, they do 
show knotty structure due to oscillations around pressure equi- 
librium with the ambient medium or to travelling perturbations. 
In both cases and also if perturbations are induced by irregu- 
larities in th e surrounding medi um (typically the backflow or 
cocoon, e.g. jMizuta et al. . 2010h. t hese could couple to KH in- 
stabilities (e.g., Agudo et alll200lh . Despite this, the FRII jets 
seem to be collimated up to distances of hundreds of kilopar- 
secs. To our knowledge, there are only two proposed mecha- 
nisms that could keep jets collimated: (i) the action of mag- 
netic fields, a mechanism which does not work if the jets are 
not magnetically dominate d, as is thought to be the case at 
kiloparsec scales (see, e.g., ISikora et al.L 120051) : or (ii) by in- 



ertial confinement, which seems not to be possible within the 
low-densi ty cavities wi t h decreasing pressure (see, e.g . simu- 
lations bv lScheck et all l2002t IPerucho & Martil 120071) where 
large-scale jets propagate. We propose here a third mechanism, 
purely hydrodynamical, that prevents the decoUimation of rel- 
ativistic jets, based on the growth of resonant modes. 

According to our results (Papers I and II and this work), this 
mechanism is more effective in jets with larger ffow Lorentz 
factors, but also depends on the initial width of the shear layer 
and on the temperature of the jet. Solutions to the linear per- 



turbation problem in a sheared, relativistic jet show that, if the 
jet is separated from the ambient by a pure contact disconti- 
nuity (vortex sheet) or by a very thick shear layer (> 1 Rjii 
the resonances do not develop, thus leading to unstable con- 
figurations in which the jets can be disrupted by the growth of 
instabilities. The same solutions, along with numerical simula- 
tions, show that the resonances grow faster in colder jets, at a 
constant Lorentz factor. In addition, it was shown in previous 
works that lower density ratios between the jet and the ambient 
medium may also lead to the disruption of the jet by the growth 
of instabilities. 

We know at lea st one case in which a jet classified as 



an FRII (0836+710. iHummel et all 119921) for its luminosity, 



seems to be disrupted by the growth of instabilities on the kilo- 
parsec s cales. This could be assign ed to the growth of insta- 
bilities (IPerucho & Lobanovi I2007D . Perhaps in this case, the 
conditions in the jet (e.g., thickness of the shear layer or phys- 
ical parameters in the jet) are such that resonant modes are not 
excited. 

One observational feature of the growth of these modes 
would be the edge brightening of jets, which should mainly 
occur in FRII jets. Possible evidence of hot shear layers is 
given by Swain, Bridle & Baum (1998), who show that the jets 
in the FRII radio galaxy 3C353 have some transversal struc- 
ture with higher radio-emission towards the boundaries of the 
jet. These hot shear layers could be generated by the growth 
of resonant modes as shown by our simulations. In addition. 



Kadler et al.l (120081) suggest the growth of resonant modes as 



the mechanism that provides the long-term stability of the jet 
in the FRII radio galaxy 3C111. A possible test for this hy- 
pothesis would be to check the transversal structure in this and 
other FRII jets seeking enhanced emission at the jet bound- 
aries. From the point of view of particle acceleration, the hot 
and turbulent (down to numerical resolution) mixing layers 
produced by the resonant modes may be sites of efficient ac- 
celeration of non-thermal par ticles (e.g.. lStawarz & Ostrowski , 



20021: lAlov & Mimicall2008l) . 



5. Conclusions 

In this work, we presented the first solutions of the linear prob- 
lem for 3D sheared relativistic flows. Following this result, we 
presented numerical simulations of the growth of KH instabil- 
ity modes from the linear to the nonlinear regimes, for three 
selected relativistic jet models. These numerical simulations 
were performed with the new high-resolution shock-capturing 
3D RHD code Ratpenat, which was parallelized using a hy- 
brid scheme with OMP and MPI processes. The simulations of 
jets in the typical stability regions of the space of parameters 
(Paper I) have confirmed previous results obtained in the 2D 
approximation. We proved that under certain conditions, the 
appearance of the shear-layer resonant modes of KH instability 
has important consequences for the stability of relativistic jets. 
In particular, the implications for the collimation of FRII jets 
were discussed. We also show, as a by-product of our simula- 



' Taking the jet radius as the distance where the velocity is half that 
in the spine, typically == 0.5 c 
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tions, that the development of different modes with a variety of 
wavelengths may induce rich transversal structure in jets, with 
shorter modes showing up in the central regions of the jet and 
longer ones dominating the structure in the m ixing layer. In 



this re spect, we confirm the results obtained in iPerucho et al 
(120061 

Future work along this line will consider the extension of 
our study in both the linear and nonlinear regime by means of 
simulations with higher resolution to a wider sample of jet pa- 
rameters, including magnetic fields and a more realistic equa- 
tion of state. 
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